home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_4.8 / thinw / thinw.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  24.6 KB  |  823 lines
  1. /* 
  2.  * thinw.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /* THINW:       program performs kxk thinning on image and yields
  10.  *            thinned image whose values are equal to the original
  11.  *              line width around that point
  12.  *              usage: thinw inimg outimg [-k MASK_SIZE] [-d] [-I] [-L]
  13.  */
  14.  
  15. #include <stdio.h>
  16. #include <string.h>
  17. #include <stdlib.h>
  18. #include <math.h>
  19. #include <images.h>
  20. #include <tiffimage.h>          /* tiff file format info */
  21. extern void print_sos_lic ();
  22.  
  23. #define DFLTMAXK 3              /* dflt maximum thinning kernel size */
  24. #define MAXMAXK 21              /* max of the maximum thinning kernel size */
  25. #define MAXITER 20              /* maximum number of iterations */
  26. #define MAXDISPLAY 40           /* maximum length of line for display */
  27.  
  28. #define OFF 0                   /* initial and final value of OFF pixels */
  29. #define ON 255                  /* initial value of ON pixels */
  30. #define ERASEBASE 51            /* OFF (or formerly erased): 0-ERASEBASE */
  31. #define ANCHORBASE 101          /* ERASED (this iteration): ERASEBASE-ANCHORBASE */
  32. #define PONBASE 152             /* ANCHOR (this iteration): ANCHORBASE-PONBASE */
  33.    /* PON (final width value): PONBASE-ON */
  34.  
  35.  
  36. long ksize (unsigned char **, struct point, long, long, long);
  37. long getring (unsigned char **, struct point, long, long, long, unsigned char *);
  38. long sqron (unsigned char **, long, long, long);
  39. long thinring (unsigned char *, long, long *);
  40. int width (unsigned char *k, long);
  41. long chkconnect (unsigned char *, long);
  42. int anchor (unsigned char *, long);
  43. int erasesqr (unsigned char **, long, long, unsigned char *, long, int, long *);
  44. long usage (short);
  45. long input (int, char **, long *, long *, long *);
  46.  
  47. main (argc, argv)
  48.      int argc;
  49.      char *argv[];
  50. {
  51.   Image *imgIO;                 /* input/output image pointer */
  52.   long maxK,                    /* max. sidelength of thinning kernel */
  53.     k,                          /* sidelength of current thinning kernel */
  54.     kM1,                        /* k - 1 */
  55.     change[MAXMAXK],            /* no. erasures for each mask size */
  56.     nIter,                      /* no. iterations */
  57.     x, y,                       /* image coordinates */
  58.     nChange,                    /* no. thinning operations on iteration */
  59.     dFlag,                      /* display results after each iter if =1 */
  60.     invertFlag;                 /* invert input image before processing */
  61.   long nONs,                    /* total ONs in original image */
  62.     nErased;                    /* cumulative no. ERASED in image */
  63.   char c;
  64.   long temp;
  65.   long permOn;                  /* permanent ON flag set to 1 if so */
  66.   unsigned char **image;        /* input/output image */
  67.   struct point imgSize;         /* image size */
  68.   unsigned char *ring;          /* ring of pixels on perimeter of kxk sqr */
  69.  
  70. /* user input */
  71.   if (input (argc, argv, &maxK, &dFlag, &invertFlag) < 0)
  72.     return (-1);
  73.  
  74. /* read input file */
  75.   imgIO = ImageIn (argv[1]);
  76.   image = imgIO->img;
  77.   imgSize.x = ImageGetWidth (imgIO);
  78.   imgSize.y = ImageGetHeight (imgIO);
  79.   if (imgSize.y > MAXDISPLAY)
  80.     dFlag = 0;
  81.  
  82. /* invert image */
  83.   if (invertFlag) {
  84.     for (y = 0; y < imgSize.y; y++)
  85.       for (x = 0; x < imgSize.x; x++)
  86.       image[y][x] = 255 - image[y][x];
  87.   }
  88.  
  89.   if ((ring = (unsigned char *) malloc (4 * (maxK - 1))) == NULL) {
  90.     printf ("not enough memory -- sorry");
  91.     return (-1);
  92.   }
  93.  
  94. /* zero image borders */
  95.   for (y = 0; y < imgSize.y; y++)
  96.     image[y][0] = image[y][imgSize.x - 1] = OFF;
  97.   for (x = 0; x < imgSize.x; x++)
  98.     image[0][x] = image[imgSize.y - 1][x] = OFF;
  99.  
  100.   for (k = 0; k < MAXMAXK; k++)
  101.     change[k] = 0;
  102.  
  103. /* iteratively convolve through image until thinned */
  104.  
  105.   nONs = nErased = 0;
  106.   for (nIter = 0, nChange = 1; (nChange > 0) && nIter <= MAXITER; nIter++) {
  107.     nChange = 0;
  108.     for (y = 1; y < imgSize.y - 1; y++) {
  109.       for (x = 1; x < imgSize.x - 1; x++) {
  110.         permOn = 0;
  111.         k = ksize (image, imgSize, x, y, maxK);
  112.         if (nIter == 0 && k >= 3)
  113.           nONs++;
  114.         kM1 = (k > 3) ? k - 1 : 3;
  115.         while (k >= kM1) {
  116.           if (sqron (image, x, y, k) == 0)
  117.             break;
  118.           if (getring (image, imgSize, x, y, k, ring) == 1)
  119.             break;
  120.           if (thinring (ring, k, &permOn) == 1) {
  121.             if (chkconnect (ring, k) == 1) {
  122.               nChange++;
  123.               (change[k])++;
  124.               erasesqr (image, x, y, ring, k, anchor (ring, k), &nErased);
  125.               break;
  126.             }
  127.           }
  128.           --k;
  129.         }
  130.         if (permOn != 0)
  131.           image[y][x] = (unsigned char) permOn;
  132.       }
  133.     }
  134.     printf ("%d interations: nChange = %7d\n", nIter + 1, nChange);
  135.     printf ("      3: %d,   4: %d,   5: %d,   6: %d,   7: %d\n\n",
  136.             change[3], change[4], change[5], change[6], change[7]);
  137.     for (y = 0; y < imgSize.y; y++) {
  138.       for (x = 0; x < imgSize.x; x++) {
  139.         if (image[y][x] >= ERASEBASE && image[y][x] < ANCHORBASE)
  140.           image[y][x] = image[y][x] - ERASEBASE;
  141.         else if (image[y][x] >= ANCHORBASE && image[y][x] < PONBASE)
  142.           image[y][x] = image[y][x] - ANCHORBASE;
  143.       }
  144.     }
  145.     if (dFlag == 1) {
  146.       for (y = 0; y < imgSize.y; y++) {
  147.         for (x = 0; x < imgSize.x; x++) {
  148.           if (image[y][x] == OFF)
  149.             printf ("  ");
  150.           else if (image[y][x] == ON)
  151.             printf ("X ");
  152.           else if (image[y][x] >= PONBASE) {
  153.             temp = (long) image[y][x] - PONBASE;
  154.             if (temp < 10)
  155.               printf ("%1d ", temp);
  156.             else
  157.               printf ("%1d", temp);
  158.           }
  159.           else {
  160.             temp = (long) image[y][x];
  161.             if (temp < 10)
  162.               printf ("%1d ", temp);
  163.             else
  164.               printf ("%1d", temp);
  165.           }
  166.         }
  167.         printf ("\n");
  168.       }
  169.       printf ("Enter <CR> for next iteration.\n");
  170.       scanf ("%c", &c);
  171.     }
  172.   }
  173.   if (dFlag == 1) {
  174.     for (y = 0; y < imgSize.y; y++) {
  175.       for (x = 0; x < imgSize.x; x++) {
  176.         if (image[y][x] == OFF)
  177.           printf ("  ");
  178.         else if (image[y][x] == ON)
  179.           printf ("X ");
  180.         else if (image[y][x] >= PONBASE) {
  181.           temp = (long) image[y][x] - PONBASE;
  182.           if (temp < 10)
  183.             printf ("%1d ", temp);
  184.           else
  185.             printf ("%1d", temp);
  186.         }
  187.         else {
  188.           printf ("- ");
  189.         }
  190.       }
  191.       printf ("\n");
  192.     }
  193.   }
  194.   if (nONs == 0)
  195.     printf ("empty image\n");
  196.   else
  197.     printf ("nONs = %d, nErased = %d (%d%%)\n", nONs, nErased,
  198.             (nErased * 100) / nONs);
  199.   if (nIter >= MAXITER && nChange != 0)
  200.     printf ("Nuts -- maximum iterations reached\n");
  201.   for (y = 1; y < imgSize.y - 1; y++)
  202.     for (x = 1; x < imgSize.x - 1; x++)
  203.       image[y][x] = (image[y][x] < PONBASE) ? OFF : ON;
  204.  
  205. /* un-invert image */
  206.   if (invertFlag) {
  207.     for (y = 0; y < imgSize.y; y++)
  208.       for (x = 0; x < imgSize.x; x++)
  209.       image[y][x] = 255 - image[y][x];
  210.   }
  211.  
  212.   ImageOut (argv[2], imgIO);
  213.  
  214.   return (0);
  215. }
  216.  
  217.  
  218. /* KSIZE:       function determines k, where kxk is largest square
  219.  *            around (x,y) which contains all ON or ERASED
  220.  *                      usage: k = ksize (image, imgSize, x, y, maxK)
  221.  */
  222.  
  223. long
  224. ksize (image, imgSize, x, y, maxK)
  225.      unsigned char **image;     /* input/output image */
  226.      struct point imgSize;      /* image size */
  227.      long x, y,                 /* image coordinates */
  228.        maxK;                    /* maximum k value */
  229. {
  230.   long k,                       /* mask size */
  231.     upHalf, downHalf,           /* half of mask below and above center */
  232.     xStart, xEnd,               /* x- start and end of square */
  233.     yStart, yEnd,               /* y- start and end of square */
  234.     xMask, yMask;               /* x,y mask coordinates */
  235.  
  236.   if (image[y][x] < ERASEBASE)
  237.     return (0);
  238.  
  239.   for (k = 4; k <= maxK; k++) {
  240.     if (k % 2 == 1)
  241.       downHalf = upHalf = (k - 3) / 2;
  242.     else {
  243.       upHalf = (k - 2) / 2;
  244.       downHalf = (k - 4) / 2;
  245.     }
  246.     xStart = x - downHalf;
  247.     xEnd = x + upHalf;
  248.     yStart = y - downHalf;
  249.     yEnd = y + upHalf;
  250.     for (yMask = yStart; yMask <= yEnd; yMask++)
  251.       for (xMask = xStart; xMask <= xEnd; xMask++)
  252.         if (image[yMask][xMask] < ERASEBASE)
  253.           return (k - 1);
  254.   }
  255.   return (maxK);
  256. }
  257.  
  258.  
  259.  
  260. /* GETRING:     function gets ring of pixels on perimeter of k-size square
  261.  *                usage: allOnes = getring (image, imgSize, x, y, k, ring)
  262.  *              If ring includes only ON pixels, function returns 1,
  263.  *              otherwise 0.
  264.  */
  265.  
  266. long
  267. getring (image, imgSize, x, y, k, ring)
  268.      unsigned char **image;     /* input/output image */
  269.      struct point imgSize;      /* image size */
  270.      long x, y,                 /* image coordinages */
  271.        k;                       /* square sidelength of ring */
  272.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  273. {
  274.   long upHalf, downHalf,        /* half of mask below and above center */
  275.     xStart, xEnd,               /* x- start and end of square */
  276.     yStart, yEnd,               /* y- start and end of square */
  277.     i, allOnes;                 /* =1 if all ones in ring; 0 otherwise */
  278.  
  279.   if (k % 2 == 1)
  280.     downHalf = upHalf = (k - 1) / 2;
  281.   else {
  282.     upHalf = k / 2;
  283.     downHalf = (k - 2) / 2;
  284.   }
  285.   xStart = x - downHalf;
  286.   xEnd = x + upHalf;
  287.   yStart = y - downHalf;
  288.   yEnd = y + upHalf;
  289.  
  290.   allOnes = 1;
  291.   i = 0;
  292.   for (x = xStart, y = yStart; x <= xEnd; x++) {
  293.     if (image[y][x] < ERASEBASE)
  294.       allOnes = 0;
  295.     ring[i++] = image[y][x];
  296.   }
  297.   for (y = yStart + 1, x = xEnd; y <= yEnd; y++) {
  298.     if (image[y][x] < ERASEBASE)
  299.       allOnes = 0;
  300.     ring[i++] = image[y][x];
  301.   }
  302.   for (x = xEnd - 1, y = yEnd; x >= xStart; --x) {
  303.     if (image[y][x] < ERASEBASE)
  304.       allOnes = 0;
  305.     ring[i++] = image[y][x];
  306.   }
  307.   for (y = yEnd - 1, x = xStart; y > yStart; --y) {
  308.     if (image[y][x] < ERASEBASE)
  309.       allOnes = 0;
  310.     ring[i++] = image[y][x];
  311.   }
  312.  
  313. /* if this square is already at border, cannot go to larger square */
  314.   if (xStart <= 0 || yStart <= 0
  315.       || xEnd >= (imgSize.x - 1) || yEnd >= (imgSize.y - 1))
  316.     return (0);
  317.  
  318.   return (allOnes);
  319. }
  320.  
  321.  
  322. /* SQRON:       function tests pixels in kxk square are already erased or 
  323.  *            for 3x3 if they can never be erased
  324.  *                      usage: flag = sqron (image, x, y, k)
  325.  *              flag = 0 if cannot erase any pixels in square
  326.  *                   = 1 if at least one pixel is still ON
  327.  */
  328.  
  329. long
  330. sqron (image, x, y, k)
  331.      unsigned char **image;     /* input/output image */
  332.      long x, y,                 /* image coordinages */
  333.        k;                       /* square sidelength of ring */
  334. {
  335.   long upHalf, downHalf,        /* half of mask below and above center */
  336.     yStart, yEnd,               /* bounds of center erase area */
  337.     xStart, xEnd;
  338.  
  339. /* check for 3x3 */
  340.   if (k == 3) {
  341.     if (image[y][x] == ON)
  342.       return (1);
  343.     else
  344.       return (0);
  345.   }
  346.  
  347. /* check center square */
  348.   if (k % 2 == 1)
  349.     downHalf = upHalf = (k - 3) / 2;
  350.   else {
  351.     upHalf = (k - 2) / 2;
  352.     downHalf = (k - 4) / 2;
  353.   }
  354.   xStart = x - downHalf;
  355.   xEnd = x + upHalf;
  356.   yStart = y - downHalf;
  357.   yEnd = y + upHalf;
  358.   for (y = yStart; y <= yEnd; y++)
  359.     for (x = xStart; x <= xEnd; x++)
  360.       if (image[y][x] >= PONBASE)
  361.         return (1);
  362.   return (0);
  363. }
  364.  
  365.  
  366.  
  367. /* THINRING:    function makes decision to thin or not based on CNUM
  368.  *            and FNUM in perimeter ring
  369.  *                      usage: flag = thinring (ring, k, &permOn)
  370.  *              Flag = 1 if thinning conditions met, 0 otherwise.
  371.  */
  372.  
  373. long
  374. thinring (ring, k, permOn)
  375.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  376.      long k;                    /* square sidelength of ring */
  377.      long *permOn;              /* permanent-ON width value */
  378. {
  379.   long nRing,                   /* no. pixels in ring */
  380.     cNum,                       /* connectivity number */
  381.     fNum,                       /* no. 1s on ring */
  382.     i, n, m;
  383.   long nOff,                    /* current run of OFF in ring */
  384.     phi0,                       /* maximum run of OFF in ring */
  385.     nFirstRun;                  /* no  OFF from ring[0] */
  386.   long lower, upper;            /* adjacent ring elements for cNum calc */
  387.  
  388.   nRing = 4 * k - 4;
  389.  
  390. /* calculate FNUM */
  391.   for (i = 0, fNum = 0; i < nRing; i++)
  392.     if (ring[i] >= ERASEBASE)
  393.       fNum++;
  394.  
  395. /* calculate 4-connected run of 0s */
  396.   nOff = (ring[0] < ERASEBASE) ? 1 : 0;
  397.   nFirstRun = (nOff == 1) ? 0 : -1;
  398.   for (i = 1, phi0 = 0; i < nRing; i++) {
  399.     if (ring[i] < ERASEBASE)
  400.       nOff++;
  401.     else {
  402.       if (nOff > 0) {
  403.         phi0 = (nOff > phi0) ? nOff : phi0;
  404.         if (nFirstRun == 0)
  405.           nFirstRun = nOff;
  406.         nOff = 0;
  407.       }
  408.     }
  409.   }
  410.   if (nOff > 0) {
  411.     if (nFirstRun > 0)
  412.       nOff += nFirstRun;
  413.     phi0 = (nOff > phi0) ? nOff : phi0;
  414.   }
  415.  
  416. /* CNUM */
  417. /* CNUM skipping corners */
  418.   for (i = 2, cNum = 0; i < nRing; i++) {
  419.     lower = (long) ring[i - 1];
  420.     if ((i % (k - 1)) == 0)
  421.       i++;                      /* skip the corner pixels */
  422.     upper = (long) ring[i];
  423.     if (upper >= ERASEBASE && lower < ERASEBASE)
  424.       cNum++;
  425.   }
  426.   if (ring[1] >= ERASEBASE && ring[nRing - 1] < ERASEBASE)
  427.     cNum++;
  428. /* CNUM at corners */
  429.   for (n = 1; n < 4; n++) {
  430.     m = n * (k - 1);
  431.     if (ring[m] >= ERASEBASE)
  432.       if (ring[m - 1] < ERASEBASE && ring[m + 1] < ERASEBASE)
  433.         cNum++;
  434.   }
  435.   if (ring[0] >= ERASEBASE && ring[1] < ERASEBASE
  436.       && ring[nRing - 1] < ERASEBASE)
  437.     cNum++;
  438.  
  439. /* to thin or not to thin */
  440.   if (cNum == 1)
  441.     if (phi0 > (k - 2) && fNum > (k - 2))
  442.       return (1);
  443.  
  444. /* for 3x3, set flag for perm. ON pixel if connection, end, or cross pt */
  445.   if (k == 3)
  446.     if (cNum > 1 || fNum <= 1 || (cNum == 0 && fNum == 4))
  447.       *permOn = width (ring, nRing) + PONBASE;
  448.  
  449.   return (0);
  450. }
  451.  
  452.  
  453. /* WIDTH:       function calculates average width from boundary of
  454.  *            permanent-ON pixel
  455.  *                      usage: w = width (ring, nRing);
  456.  *              Width calculation is made by first finding minima of
  457.  *              each ERASED run. Then the minima are averaged and
  458.  *              multiplied by 2 to obtain the width estimate.
  459.  */
  460.  
  461. #define MAXRUNSERASED 4         /* max. no. runs of erased pix.s (for cross) */
  462.  
  463. int
  464. width (ring, nRing)
  465.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  466.      long nRing;                /* no. pixels in ring */
  467. {
  468.   long min[MAXRUNSERASED],      /* minimum distance for each erased run */
  469.     nMin,                       /* no. runs of ERASED */
  470.     width,                      /* estimated width measurement */
  471.     fNum,                       /* no. ONs in nbrhood */
  472.     i, j;
  473.  
  474. /* initialize minimas */
  475.   for (i = 0; i < MAXRUNSERASED; i++)
  476.     min[i] = 255;
  477.  
  478. /* calculate minima for each run of OFFs in nbrhood */
  479.   for (i = 0, nMin = 0, fNum = 0; i < nRing; i++) {
  480.     if (ring[i] < ERASEBASE) {
  481.       for (j = i; ring[j] < ERASEBASE && j < nRing; j++) {
  482.         if (ring[j] < min[nMin])
  483.           min[nMin] = (long) ring[j];
  484.         else if (ring[j] > PONBASE)
  485.           fNum++;
  486.       }
  487.       nMin++;
  488.       i = j;
  489.     }
  490.   }
  491.   if (ring[0] < ERASEBASE && ring[nRing - 1] < ERASEBASE && fNum > 0) {
  492.     min[0] = (min[0] < min[nMin - 1]) ? min[0] : min[nMin - 1];
  493.     --nMin;
  494.   }
  495.  
  496. /* calculate average width measurement */
  497.   switch (nMin) {
  498.   case 1:
  499.     width = 2 * min[0] + 1;
  500.     break;
  501.   case 2:
  502.     width = min[0] + min[1] + 1;
  503.     break;
  504.   case 3:
  505.   case 4:
  506.     for (i = 0, width = 0; i < nMin; i++)
  507.       width += min[i];
  508.     width = (long) ((2.0 / nMin) * width + 0.5) + 1;
  509.     break;
  510.   }
  511.   if (width == 0)
  512.     width = 1;
  513.   return (width);
  514. }
  515.  
  516.  
  517.  
  518. /* CHKCONNECT:  function checks connectivity to see if current erasures
  519.  *            combined with past erasures will destroy connectivity
  520.  *                      usage: flag = chkconnect (ring, k)
  521.  *              Function returns flag = 0 if connectivity destroyed,
  522.  *              flag = 1 if connectivity retained.
  523.  */
  524.  
  525. /* if corner or its nbr is ON, then that is value of side, otherwise OFF */
  526. #define CORNER(ISIDE,ICORNER,ICORNERSIDE,ICORNEROTHER) \
  527.         if (ring[ICORNER] >= PONBASE \
  528.              || ring[ICORNEROTHER] >= ERASEBASE) \
  529.               side[ISIDE] = ON; \
  530.         else side[ISIDE] = OFF
  531.  
  532. /* for a NW corner, check if corner is anchor point, and if so set to ON */
  533. #define CORNERNW(ISIDE,ICORNER,ICORNERSIDE,ICORNEROTHER) \
  534.         if (ring[ICORNER] >= ANCHORBASE \
  535.             || ring[ICORNEROTHER] >= ERASEBASE) \
  536.               side[ISIDE] = ON; \
  537.         else side[ISIDE] = OFF
  538.  
  539. long
  540. chkconnect (ring, k)
  541.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  542.      long k;                    /* square sidelength of ring */
  543. {
  544.   unsigned char *side;          /* array of 8-connected side pixels */
  545.   long nRing,                   /* no. pixels in ring */
  546.     nSide,                      /* no. pixels along 8-connected side */
  547.     iSide,                      /* side array index */
  548.     anON,                       /* an ON run along side */
  549.     ONandERASE,                 /* run of ON and n ERASEs */
  550.     N1,                         /* no. ones in nbrhood */
  551.     i;
  552.  
  553.   nRing = 4 * k - 4;
  554.  
  555. /* calculate N1 */
  556.   for (i = 0, N1 = 0; i < nRing; i++)
  557.     if (ring[i] >= PONBASE)
  558.       N1++;
  559.   if (N1 == 0)
  560.     return (0);
  561.  
  562.   nSide = k;
  563.   if ((side = (unsigned char *) malloc (nSide)) == NULL) {
  564.     printf ("not enough memory -- sorry");
  565.     return (-1);
  566.   }
  567.  
  568. /* check connectivity of west side */
  569.   i = 3 * (k - 1);
  570.   CORNER (0, i, i + 1, i - 1);
  571.   for (i = i + 1, iSide = 1; i < nRing; i++, iSide++)
  572.     side[iSide] = ring[i];
  573.   CORNERNW (nSide - 1, 0, nRing - 1, 1);
  574.  
  575.   for (i = 0, anON = 0, ONandERASE = 0; i < nSide; i++) {
  576.     if (side[i] >= PONBASE) {
  577.       if (ONandERASE == 1) {
  578.         free (side);
  579.         return (0);
  580.       }
  581.       else
  582.         anON = 1;
  583.     }
  584.     else if ((side[i] >= ERASEBASE && side[i] < PONBASE)
  585.              && anON == 1)
  586.       ONandERASE = 1;
  587.     else if (side[i] < ERASEBASE)
  588.       anON = ONandERASE = 0;    /* off */
  589.   }
  590.  
  591. /* check connectivity of north side */
  592.   i = k - 1;
  593.   CORNER (0, i, i - 1, i + 1);
  594.   for (i = i - 1, iSide = 1; i >= 0; --i, iSide++)
  595.     side[iSide] = ring[i];
  596.   CORNERNW (nSide - 1, 0, 1, nRing - 1);
  597.  
  598.   for (i = 0, anON = 0, ONandERASE = 0; i < nSide; i++) {
  599.     if (side[i] >= PONBASE) {
  600.       if (ONandERASE == 1) {
  601.         free (side);
  602.         return (0);
  603.       }
  604.       else
  605.         anON = 1;
  606.     }
  607.     else if ((side[i] >= ERASEBASE && side[i] < PONBASE)
  608.              && anON == 1)
  609.       ONandERASE = 1;
  610.     else if (side[i] < ERASEBASE)
  611.       anON = ONandERASE = 0;    /* off */
  612.   }
  613.  
  614. /* check connectivity of east side */
  615.   i = k - 1;
  616.   CORNER (0, i, i + 1, i - 1);
  617.   for (i = i + 1, iSide = 1; iSide < nSide - 1; i++, iSide++)
  618.     side[iSide] = ring[i];
  619.   CORNER (nSide - 1, 2 * (k - 1), 2 * (k - 1) - 1, 2 * (k - 1) + 1);
  620.  
  621.   for (i = 0, anON = 0, ONandERASE = 0; i < nSide; i++) {
  622.     if (side[i] >= PONBASE) {
  623.       if (ONandERASE == 1) {
  624.         free (side);
  625.         return (0);
  626.       }
  627.       else
  628.         anON = 1;
  629.     }
  630.     else if ((side[i] >= ERASEBASE && side[i] < PONBASE)
  631.              && anON == 1)
  632.       ONandERASE = 1;
  633.     else if (side[i] < ERASEBASE)
  634.       anON = ONandERASE = 0;    /* off */
  635.   }
  636.   return (1);
  637. }
  638.  
  639.  
  640. /* ANCHOR:      function returns 1 if core is exposed on NW side
  641.  *                    usage: anchor (ring, k)
  642.  */
  643.  
  644. int
  645. anchor (ring, k)
  646.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  647.      long k;                    /* square sidelength of ring */
  648. {
  649.   long nRing,                   /* no. pixels in ring */
  650.     i;
  651.  
  652.   nRing = 4 * k - 4;
  653.   for (i = 0; i < k; i++)
  654.     if (ring[i] >= ERASEBASE)
  655.       return (0);
  656.   for (i = 3 * (k - 1) + 1; i < nRing; i++)
  657.     if (ring[i] >= ERASEBASE)
  658.       return (0);
  659.   return (1);
  660. }
  661.  
  662.  
  663.  
  664. /* ERASESQR:    function erases square contained within square perimeter
  665.  *            to minimum width values from the edge
  666.  *                      usage: erasesqr (image, x, y, ring, k, anchor, 
  667.  *                                                              &nErased)
  668.  *              If the core is an anchor, then the pixels are erased
  669.  *              to ERASED + 1, otherwise they are erased to ERASED.
  670.  *              For kxk > 3x3 erasure, PON pixels (permanent ON for 3x3)
  671.  *              are erased; and if an anchor point (ERASED + 1) can be
  672.  *              erased to a non-anchor point by a larger k (I don't know
  673.  *              if this can actually happen... I'm pretty sure it cannot
  674.  *              and have found on a fingerprint it does not, but have not
  675.  *              thought it out yet.) then it is erased to ERASED.
  676.  */
  677.    /* want to erase if ERASED +1, not just ON */
  678. int
  679. erasesqr (image, x, y, ring, k, anchor, nErased)
  680.      unsigned char **image;     /* input/output image */
  681.      long x, y,                 /* image coordinages */
  682.        k;                       /* square sidelength of ring */
  683.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  684.      int anchor;                /* 1 if core is NW pt of NW-SE diag; else 0 */
  685.      long *nErased;             /* no. of erased */
  686. {
  687.   long nRing,                   /* no. pixels in ring */
  688.     dist;                       /* minimum dist in ring */
  689.   long upHalf, downHalf,        /* half of mask below and above center */
  690.     yStart, yEnd,               /* bounds of center erase area */
  691.     xStart, xEnd, i;
  692. /* find minimum distance value in neighborhood */
  693.   nRing = 4 * k - 4;
  694.   dist = 256;
  695.   for (i = 0; i < nRing; i++)
  696.     if ((long) ring[i] < dist)
  697.       dist = (long) ring[i];
  698.   dist += k - 2;
  699.  
  700. /* erase for 3x3 */
  701.   if (k == 3) {
  702.     if (image[y][x] == ON) {
  703.       (*nErased)++;
  704.       image[y][x] = (anchor == 0) ? (unsigned char) (ERASEBASE + dist) :
  705.         (unsigned char) (ANCHORBASE + dist);
  706.     }
  707.   }
  708. /* erase for kxk > 3x3 */
  709.   else {
  710.     if (k % 2 == 1)
  711.       downHalf = upHalf = (k - 3) / 2;
  712.     else {
  713.       upHalf = (k - 2) / 2;
  714.       downHalf = (k - 4) / 2;
  715.     }
  716.     xStart = x - downHalf;
  717.     xEnd = x + upHalf;
  718.     yStart = y - downHalf;
  719.     yEnd = y + upHalf;
  720.     for (y = yStart; y <= yEnd; y++) {
  721.       for (x = xStart; x <= xEnd; x++) {
  722.         if (image[y][x] >= ANCHORBASE) {
  723.           if (image[y][x] >= PONBASE)
  724.             (*nErased)++;
  725.           image[y][x] = (anchor == 0) ? (unsigned char) (ERASEBASE + dist) :
  726.             (unsigned char) (ANCHORBASE + dist);
  727.         }
  728.       }
  729.     }
  730.   }
  731.   return (0);
  732. }
  733.  
  734.  
  735.  
  736.  
  737.  
  738. /* USAGE:       function gives instructions on usage of program
  739.  *                    usage: usage (flag)
  740.  *              When flag is 1, the long message is given, 0 gives short.
  741.  */
  742.  
  743. long
  744. usage (flag)
  745.      short flag;                /* flag =1 for long message; =0 for short message */
  746. {
  747.  
  748. /* print short usage message or long */
  749.   printf ("USAGE: thinw inimg outimg [-k MASK_SIZE] [-d] [-I] [-L]\n");
  750.   if (flag == 0)
  751.     return (-1);
  752.  
  753.   printf ("\nthinw performs iterative thinning of binary objects in input\n");
  754.   printf ("image to produce skeleton image with values OFF (0) ON (255);\n");
  755.   printf ("skeleton pixel values are widths of the original image lines.\n\n");
  756.   printf ("ARGUMENTS:\n");
  757.   printf ("    inimg: input image filename (TIF)\n");
  758.   printf ("   outimg: output image filename (TIF)\n\n");
  759.   printf ("OPTIONS:\n");
  760.   printf ("  -k MASK_SIZE: window size for kxk mask;(k>=3, default = %d\n", DFLTMAXK);
  761.   printf ("            -d: display results of each iteration\n");
  762.   printf ("                this only displays for images <= 40x40 to\n");
  763.   printf ("                fit on screen.\n");
  764.   printf ("            -I: invert input image before processing\n");
  765.   printf ("            -L: print Software License for this module\n");
  766.  
  767.   printf ("                NOTE: the image output values are the following:\n");
  768.   printf ("                -- thinned lines -- width value plus offset of %3d\n", PONBASE);
  769.   printf ("                -- off -- %3d\n", OFF);
  770.   printf ("                -- on -- %3d\n", ON);
  771.  
  772.   return (-1);
  773. }
  774.  
  775.  
  776. /* INPUT:       function reads input parameters
  777.  *                  usage: input (argc, argv, &maxK, &dFlag)
  778.  */
  779.  
  780. #define USAGE_EXIT(VALUE) {usage (VALUE); return (-1);}
  781.  
  782. long
  783. input (argc, argv, maxK, dFlag, invertFlag)
  784.      int argc;
  785.      char *argv[];
  786.      long *maxK,                /* max. sidelength of thinning kernel */
  787.       *dFlag,                   /* display results after each iter if =1 */
  788.       *invertFlag;              /* invert input image before processing */
  789.  
  790. {
  791.   long n;
  792.  
  793.   if (argc == 1)
  794.     USAGE_EXIT (1);
  795.   if (argc == 2)
  796.     USAGE_EXIT (0);
  797.  
  798.   *maxK = DFLTMAXK;
  799.   *dFlag = 0;
  800.   *invertFlag = 0;
  801.  
  802.   for (n = 3; n < argc; n++) {
  803.     if (strcmp (argv[n], "-k") == 0) {
  804.       if (++n == argc || argv[n][0] == '-')
  805.         USAGE_EXIT (0);
  806.       *maxK = atol (argv[n]);
  807.     }
  808.     else if (strcmp (argv[n], "-d") == 0)
  809.       *dFlag = 1;
  810.     else if (strcmp (argv[n], "-I") == 0)
  811.       *invertFlag = 1;
  812.     else if (strcmp (argv[n], "-L") == 0) {
  813.       print_sos_lic ();
  814.       exit (0);
  815.     }
  816.     else
  817.       USAGE_EXIT (0);
  818.   }
  819.   if (*maxK < 3)
  820.     USAGE_EXIT (1);
  821.   return (0);
  822. }
  823.